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The quantum dynamics of the linearly polarized Gowdy T 3 model (compact inhomogeneous uni- 
verses admitting linearly polarized gravitational waves) is analyzed within Loop Quantum Cos- 
mology by means of an effective dynamics. The analysis, performed via analytical and numerical 
methods, proves that the behavior found in the evolution of vacuum (homogeneous) Bianchi 1 uni- 
verses is preserved qualitatively also in the presence of inhomogeneities. More precisely, the initial 
singularity is replaced by a big bounce which joins deterministically two large classical universes. 
In addition, we show that the size of the universe at the bounce is at least of the same order of 
magnitude (roughly speaking) as the size of the corresponding homogeneous universe obtained in 
the absence of gravitational waves. In particular, a precise lower bound for the ratio of these two 
sizes is found. Finally, the comparison of the amplitudes of the gravitational wave modes in the 
distant future and past shows that, statistically (i.e., for large samples of universes), the difference 
in amplitude is enhanced for nearly homogeneous universes, whereas this difference vanishes in inho- 
mogeneity dominated cases. The presented analysis constitutes the first systematic effective study 
of an inhomogeneous system within Loop Quantum Cosmology, and it proves the robustness of the 
results obtained for homogeneous cosmologies in this context. 
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I. INTRODUCTION 

Loop Quantum Gravity (LQG) [3, is one the most 
promising approaches for the quantization of gravity. Its 
development in recent years provides hopes in the pro- 
gram to overcome the limits of Classical Relativity, like 
for example in avoiding the breakdown of physics at 
spacetime singularities. In part, this task has already 
been achieved in the context of the application of LQG 
methods to the study of (simple) cosmological systems, 
known as Loop Quantum Cosmology (LQC) [3j. In the 
simplest (isotropic) models, the initial big bang singu- 
larity has been shown to be resolved dynamically [J- 
6J, being replaced by a big bounce connecting an ex- 
panding universe with a preceding (contracting) one in 
a deterministic way. This property, initially described in 
flat Friedmann-Robertson- Walker universes with a mass- 
less scalar field, was next confirmed analytically [7[ as 
well as extended to models of spherical and hyperbolic 
topology Q or with nontrivial cosmological constant Q , 
and even further to anisotropic homogeneous systems in 
vacuo pjJGIl or with matter [IMl]. 

Even within these simple settings, the quantum ef- 
fects change drastically the standard understanding of 
the early universe cosmology, in particular providing so- 
lutions to the horizon problem, conserving the (otherwise 
violated) entropy bounds [15[ , or increasing the estimate 
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on the probability of inflation as to become an almost 
certain event [la ]. It is also expected that the singular- 
ity resolution mechanism of LQC extends to black hole 
interiors [17| , extension which may allow one to cure the 
long standing problems of black hole physics (for exam- 
ple preventing potential information loss) and reveal yet 
unpredicted phenomena. 

In spite of this success, the commented mechanism to 
avoid singularities cannot be taken for granted yet since 
it needs further confirmation in more realistic situations, 
which necessarily include inhomogeneities. The few at- 
tempts to achieve this goal that have appeared in the lit- 
erature [Lsl [l9| rely on heuristic constructions, inspired 
by a polymer quantization, but are not based on any 
quantum model that can be considered rigorously de- 
fined. 

Until recently, the task of building a model out of 
a well defined quantum theory, with good understand- 
ing and control on each of the steps of its construc- 
tion, seemed out of reach, as the most promising midi- 
superspace treatments [2(| were never developed past the 
quantum kinematics and the formal introduction of the 
constraints. In the past few years, however, an alternate 
possibility emerged when a quantization scheme within 
the framework of LQC was formulated [21-23] for a class 
of spatially compact cosmological spacetimes admitting 
linearly polarized gravitational waves: the T 3 Gowdy 
universes [Hj]. The spacetimes of this class admit two 
spatial Killing fields (with space of orbits diffcomorphic 
to S 1 ). Since they are highly symmetric, a symmetry 
reduction is allowed. But, even so, they still possess lo- 
cal degrees of freedom and, furthermore, their structure 
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makes a nonperturbative analysis viable. 

In order to describe these universes quantum mechani- 
cally, a so-called hybrid quantization scheme was applied 
(see Ref. [22J for details). First, the geometry was repre- 
sented as the Fourier modes of gravitational waves living 
on a homogeneous (Bianchi I O E3, [H, HBl ) background. 
Next, the background geometry was quantized using loop 
techniques, while a standard Fock quantization was used 
to represent the gravitational wave modes [13] . 

Such a quantization prescription allows one to rigor- 
ously describe the system on a kinematical level, and 
identify from it a Hilbert space of physical states [22[ 
following Dirac's program. However, the complexity of 
the quantum constraints and the presence of an infinite 
number of degrees of freedom, as well as problems related 
to the choice of an internal time in vacuo ll], make the 
analysis of the genuine quantum dynamics extremely dif- 
ficult. One of possibilities to deal with these difficulties 
is to divide the program into two steps: 

(i) A detailed analysis of the dynamical properties of 
the system by means of the so-called effective dy- 
namics [28( | , constructed on the basis of a rigorously 
defined quantum model. 

(ii) The confirmation of the validity of the effective pre- 
dictions via a genuine quantum analysis performed 
in a computationally manageable domain. 

The treatment via effective dynamics, while allowing one 
to adopt the procedures and methods of classical mechan- 
ics, incorporates the quantum effects of the discrete ge- 
ometry to a certain extent. In the cases tested so far, this 
kind of effective dynamics has proven to reproduce the 
genuine quantum dynamics of states with a semiclassical 
behavior up to a remarkable accuracy §, 0, El, EH ■ 
On the other hand, the results of this effective approach, 
even when regarded as preliminary, make accessible and 
unveil properties of the system which are important for 
the correct formulation and refinement of the genuine 
theory. 

In this article, we focus on step (0. Starting from 
the quantum model obtained with the hybrid quantiza- 
tion, we construct the effective, classical dynamical sys- 
tem which approximates the quantum dynamics to first 
order (i.e. without taking into account the effect of state- 
dependent parameters on the dynamics). We then apply 
the resulting description to address the following ques- 
tions. Are the cosmological singularities dynamically re- 
solved by a big bounce mechanism in the presence of the 
(inhomogeneous) gravitational waves? And, if the answer 
is in the affirmative, does this bounce happen under con- 
ditions which are similar to those found in homogenous 
scenarios? Or do there exist substantial modifications? 
Furthermore, if the bounce persists, how does it affect 
the gravitational wave modes? In particular, we are go- 
ing to focus our attention on the problem of whether 
the energy distribution of the modes change significantly 
through the bounce process. 



The main conclusions of this study were briefly re- 
ported in Ref. [29]. In this article, we present the de- 
tails of the analysis and the methodology employed in 
it. This material, while important to validate the re- 
sults of Ref. [1^], also provides tools which are useful 
in the analysis of more complicated (and more realistic) 
inhomogeneous cosmological systems, e.g. the corrected 
vacuum Gowdy universe description of Ref. 12311 , as well 
as its extension incorporating matter fields [301 ] . The re- 
sults themselves constitute only a qualitative guideline 
to elucidate the kind of phenomena that are observed in 
the considered class of inhomogeneous systems. The fine 
details of these results must not be regarded as provid- 
ing definitive answers, in particular owing to the specific 
quantization prescription that has been adopted (see the 
discussion in Sec. Ill Bp . The system obtained with the 
corrected quantization prescription of Ref. (23|, however, 
has a sufficiently similar structure [23| and behavior [l4| 
as to expect that our results are still applicable to it, at 
least qualitatively. 

The paper is organized as follows. We start by pro- 
viding a brief introduction to the hybrid quantization of 
the model in Sec. HU Next, in Sec. [En] we introduce its 
effective description, deriving and discussing the struc- 
ture of the effective equations of motion (EOM) and their 
asymptotic behavior. In order to provide a reference for 
comparison, these equations are analytically solved in 
the absence of inhomogeneities in Sec. IIV1 thus integrat- 
ing the effective dynamics of the Bianchi I homogeneous 
background. The core of the article is presented in two 
sections. In Sec. [V] we combine analytical and numerical 
methods to discuss how the presence of inhomogeneities 
affects the appearance and properties of the bounce. On 
the other hand, Sec. IVII is devoted to numerical studies 
of the dynamical trajectories. Finally, in Sec. I VIII we 
study the effect of the bounce on the structure of the in- 
homogeneities by Monte-Carlo methods. In Sec. lVIlH we 
provide a summary of the results, as well as a discussion 
of their relevance, limitations, and possible extensions. 

II. HYBRID GOWDY MODEL 

In this section we briefly introduce the class of Gowdy 
universes of T 3 topology quantized within the LQC 
framework. We start by specifying the model at the clas- 
sical level in Subsec. lll Al and outline then the main steps 
and properties of the hybrid quantization in Subsec. lIIBI 

A. The Gowdy T 3 universe 

The general Gowdy models are spacetimes with spatial 
sections of compact topology and two spatial Killing vec- 
tor fields. In this article, we restrict our attention just to 
one model: the vacuum model with the spatial topology 
of a 3-torus (T 3 ). Besides, we will consider only the sub- 
family of universes whose content of gravitational waves 
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is linearly polarized. These restrictions imply that the 
Killing fields are axial and hypersurface orthogonal. 

Within this family of spacetimes, we can make use of a 
distinguished coordinate system (t,9,a,S), where a and 
6 are axial coordinates adapted to the Killing fields. In a 
3 + I decomposition, owing to the hypersurface orthogo- 
nality, the induced 3-metric can be chosen to be diagonal, 
and one characterizes it just by three fields correspond- 
ing to the norm of one of the Killing fields, the area of 
the isometry group orbits, and the scale factor of the 1- 
metric induced on the manifold of group orbits. Most of 
the gauge freedom is fixed by the condition that both the 
generator of conformal transformations of the considered 
1-metric as well as the area of the orbits be homogeneous 
functions. In fact this leaves out only two gauge degrees 
of freedom, namely, those corresponding to the spatial 
averages of the Hamiltonian constraint and of the difieo- 
morphism acting in the 9 direction. 

The phase space of the system splits into two sectors 
formed, respectively, by the following degrees of freedom 

a) Homogeneous: the spatial averages (zero Fourier 
modes) of the metric fields and their conjugate vari- 
ables. 

b) Inhomogeneous: all the nonzero modes of the 
only metric field £(0) whose spatial dependence has 
not been determined in the gauge-fixing procedure 
(namely, the norm of one of the Killing vector fields) , 
together with the conjugate momenta of these modes. 

The homogeneous sector describes the phase space of 
a vacuum Bianchi I model with 3-torus topology [3lj ]. 
to which the considered Gowdy model reduces when all 
the inhomogeneities vanish. In order to represent it, we 
employ the standard LQC description of the Bianchi I 
model [To| , expressing the degrees of freedom in terms of 
Ashtekar-Barbero variables. These are an SU(2) valued 
connection Af, and a densitized triad E l a . Upon fixing 
the gauge by the requirement that the 3-metric be di- 
agonal and with the introduction of a fiducial Euclidean 
metric, these variables get determined by (time varying) 
connection and triad coefficients: 
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Here, a is an internal SU(2) index, i = 9, a, 5, and no 
summation is assumed for repeated indices. These coef- 
ficients are canonically conjugate, with Poisson brackets 
{c l ,pj} = 8irG"f5j, where 7 is the Immirzi parameter. 

The degrees of freedom of the inhomogeneous sector 
are encoded in the metric field £(9) and its conjugate 
momentum P^(9), both with the zero mode excluded. 
They can be decomposed into (nonzero) Fourier modes, 
obtaining an infinite countable set of canonical pairs, 
{(£m, P™), m € Z \ {0}}. For these modes, we introduce 
creation and annihilation variables (a m ,a^J in analogy 
to those naturally associated with a massless scalar field: 
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It is easy to check that they have the Poisson brackets 
{a m ,a T -} — ~i5 m m- 

With this choice of variables, the spatial metric takes 
the form: 



ds 



2 _ \pePaPs\ 



4tt 2 



N 2 



d0 2 



-I21L dcr 2 



j?i£d<5 2 



+ e ^— + 6^^ 

pi p s 



where N is the densitized lapse function, 
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is the contribution of the nonzero modes to the field £(6), 
and the function y(9) equals 
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The only nontrivial constraints of the system -the spa- 
tial averages of the densitized Hamiltonian constraint Cq 
and of the diffeomorphism constraint in the 9 direction 
Cg (corresponding to rigid rotations in 9)— can be written 



Cg = ^ m ( a m a m ~ fl -m«-m); 
m— 1 
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In these equations, all the information about the inho- 
mogeneous degrees of freedom is contained in the terms 
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which have the form of a free field Hamiltonian and of an 
interaction Hamiltonian (which creates and annihilates 
pairs of particles). Finally, defining 0^ = CiP%, we can 
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express the classical Hamiltonian constraint Cq in a more 
compact form: 

c G = -^(Qe®s + e e e CT + e^) 

1 (2.8) 

+ — (0, + e ff ) 2 — Hf at + 32n 2 G\ Pe \Hl 
7 \Pe\ 

At this point, it is worth commenting that, if one 
adopts the same choice of lapse as in full LQG (namely 
N = 1), the zero mode that one obtains for the Hamil- 
tonian constraint is Cg/V, rather than Cq (with V be- 
ing the Bianchi I volume V := \/\pePsPa\)- In fact, the 
constraint Cg/V is indeed the one which is implemented 
quantum mechanically in the treatment applied in Refs. 
6, 10]. The counterpart of Cg is then constructed by 
introducing a scaling by V at the quantum level. 

B. Hybrid quantization 

The quantization of this Gowdy model was achieved 
recently adopting a hybrid approach [2TM23L l3ll ] , which 
is based on the hypothesis that the most relevant quan- 
tum geometry effects are those affecting the homogeneous 
sector representing the Bianchi I background on which 
the gravitational waves propagate. In this way, the zero 
modes of the geometry are quantized using the techniques 
of LQG, while a standard Fock quantization is employed 
for the inhomogeneities arising from the presence of grav- 
itational waves. 

The proposed quantum model is thus applicable in the 
regime where the "energy" of the inhomogeneities is dis- 
tributed among the modes so that, while the overall back- 
reaction of the background may feel the discreteness ef- 
fects, the energy of each individual mode is still small for 
the polymer effects on it to become important. Remark- 
ably, this hybrid quantization proves to be viable and 
consistent. The kinematical Hilbert space of the system 
is constructed as the product of the kinematical Hilbert 
space for Bianchi I (with T 3 -topology) in LQC and the 
Fock space for the inhomogeneous sector. The two sec- 
tors are coupled by the Hamiltonian constraint, so that 
the quantization is indeed nontrivial. 

1. The Bianchi I background 

Let us first describe the polymeric quantization of the 
Bianchi I geometry, i.e., of the homogeneous sector. In 
this quantization scheme, the elementary variables are 
fluxes, which are proportional to the triad coefficients Pi, 
and holonomies of connections computed along edges of 
coordinate length in each of the (gauge-fixed) di- 

rections i, where fj,j is any real number. The elements of 
these holonomies are linear combinations of the exponen- 
tials A/^j (fii) — exp(i/^iCi/2) (no Einstein summation con- 
vention is adopted). In the triad representation, each of 



these elements is represented by an eigenstate The 
completion of the vector space spanned by these states 
with respect to the discrete inner product (Atj|/4) = 5^1. 
provides the kinematical Hilbert space, %Kin,i, for each 
direction. The kinematical Hilbert space of the homoge- 
neous Sector is Simply %Kin — ®iHKm,i- 

On the basis of states \fJ>i), the action of the elementary 
operators is 

pilm) = iTrjl^fiilfii), A/^l^) = |/ij + /4). (2.9) 

Here, Ipi = v 7 Gh is the Planck length. The Hamil- 
tonian constraint is constructed by defining curvatures 
in terms of holonomies along closed loops, formed by 
straight edges in the different directions 8, a, and 6, and 
regularizing the inverse volume by means of commutators 
of holonomies with the volume operator 0, [25|, |32[ • 

To encode the discreteness of the geometry character- 
istic of LQG, in the above construction of the curvature 
operator in terms of holonomies, the area enclosed by the 
loops is assumed to coincide with the minimum nonzero 
eigenvalue A of the area operator in LQG. This prescrip- 
tion is usually called the improved dynamics approach. 
As a consequence, the coordinate length of the holonomy 
along each edge takes a value 2-7r/2i which depends on the 
considered state, because so does the area. 

In the process of arriving to a proper understanding of 
the anisotropic models in LQC, different ways to imple- 
ment this prescription have been presented in the litera- 
ture [HI, H3] • We will adopt here the proposal explored, 
e.g., in Refs. [l(| HEl , according to which 

I = Vp. (2.10) 

Although this proposal has proven to lead to physical in- 
consistencies in the case of noncompact spatial sections 
(see e.g. [34|), we will follow it here owing to two reasons: 

(i) its simple mathematical structure provides a well con- 
trollable way of combining the resulting background with 
the Fock theory which describes the inhomogeneities, and 

(ii) the proposal that corrects its drawbacks, put forward 
in Ref. [LI, is similar enough in structure [lj] as to ex- 
pect that the qualitative results obtained here will also 
hold for it. We also note that the commented problems 
of inconsistency are not present in our system, owing to 
the T 3 -topology [Hj], and that the possible discrepancies 
with the proposal of Ref. [l3[ should disappear in the 
region of large triad variables. 

Given that fii depends on the triad coefficients, the 
operator generates in fact a state dependent trans- 
formation on the basis formed by jj, a , /u,$) — 
Fortunately, the generators (in phase space) of the cor- 
responding shifts commute, so that one can relabel this 
basis using affine parameters Vj instead of the labels /itj, 
in such a way that the action of Afp i is given just by a 
constant displacement of the new label Vj. These affine 
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parameters are related to the piS as follows 

^sgnfe)H 3/2 - 
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With this relabeling, the action of the basic operators 
becomes [H [H 
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On the other hand, the regularized triad operator, rep- 
resenting the inverse of Pi, takes the form: 
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It is possible to see |22|, [31J that, with a suitable choice 
of symmetric factor ordering, the Hamiltonian constraint 
attained for the Bianchi I model with this quantization 
approach (and with the lapse chosen equal to the unity) 
leaves invariant the superselection sectors formed by the 
states which have support only on (the product over the 
directions i of) semilattices of the form 



C% ={±(e 4 + 2fc),fc = 0,l,2...}, 



(2.15) 



where Si can be any real number in the interval (0,2]. 
Note that the states with a vanishing Vj decouple from 
these superselection sectors. In particular, this allows one 
to introduce the inverse of the operator (|2.13l) in each of 
those sectors. In turn, this enables one to make a change 
of densitization in the Hamiltonian constraint [2l|, [22], HH 
which provides us, finally, with a quantum counterpart 
of the Bianchi I part of the constraint Cg, introduced in 
Eq. (|2.6b[) . This constraint is obtained by representing 
the classical quantities &i by the operators: 



2. The inhomogeneities 

We next proceed to introduce and quantize the inho- 
mogeneities of the system. This is achieved by construct- 
ing a symmetric Fock space J-, where the creation and an- 
nihilation variables {(a m , a^)} f° r au the nonzero modes 
are promoted to creation and annihilation operators act- 
ing on the vacuum. This Fock quantization is the only 
one (up to equivalence) which respects the invariance of 
the vacuum under rigid rotations in 9 and ensures that 
(after deparametrization) the field dynamics is unitarily 
implemented (36j | . 

The n-particle states |{n m }) := |..., n_ m , n m , ...) 
provide an orthonormal basis for this Fock space. Here 
n rn < oo is the occupation number of the m-th mode, and 
it is assumed that only a finite set of these occupation 
numbers is different from zero. 

Adopting normal ordering, the constraint Cg given in 
Eq. (|2.6ap . which is the generator of rigid rotations, be- 
comes the quantum operator 



C B = 



m(aj, 



(2.18) 



The kernel of this operator is a proper Fock subspace. 

In the total kinematical Hilbert space HKin®J~, we can 
follow the quantization procedure explained above and 
construct the quantum constraint representing (|2.6bl) . 
This constraint is obtained by adopting normal ordering 
in the operator representation of the inhomogeneous con- 
tributions (|2.7a[) and (|2.7bl) , and by substituting pg and 
&i by their quantum counterparts (|2.12bl) and (I2.16[) . 

Finally, the physical Hilbert space of the system can 
be found, e.g., by realizing that one can identify solu- 
tions of the Hamiltonian constraint with initial data on 
the section of minimum value of vg, namely vg = eg. The 
Hilbert structure in the space of initial data can be deter- 
mined by impo sing se lf-adjointness conditions on physi- 
cal observables 22 , 3l[ . Then, on the constructed Hilbert 



space, one can easily impose the remaining constraint Cg 
completing the quantization of the model. 
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Mi = 
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(2.16) 



(2.17) 



Note that all these operators commute among them- 
selves, and 0^ commutes with pj for i =/= j. 



III. EFFECTIVE GOWDY DYNAMICS 

In this section we introduce the effective EOM that 
result from the hybrid quantization of the Gowdy space- 
times, and solve them in certain regimes of interest. This 
effective treatment is going to be classical, in the sense 
that we are going to define an effective classical (poly- 
merized) Hamiltonian and the evolution of the relevant 
quantities is going to be provided by Hamilton- Jacobi 
equations. Even so, this Hamiltonian captures in fact 
relevant information about the modifications to General 
Relativity that are expected to mimic the (loop) quan- 
tum behavior of the system. 



6 



A. Equations of motion 



Given the quantum treatment of the system presented 
in Sec. Ill Bt we now proceed with the construction of its 
effective description. We ignore all the state-dependent 
quantities, thus building only a first-order description. 
This amounts to replacing the basic operators, that is Pi 
and Af±ji i [see Eq. (|2.12[) ]. as well as a m and a^, by their 
expectation values. This constitutes an approximation 
which should be valid for very sharply peaked states. To 
simplify the notation, we will further drop the brackets 
(•) in the expectation values, denoting for example (pi) 
by pi, and so on. 

Furthermore, since we are interesting in solutions cor- 
responding to macroscopic universes, we focus our atten- 
tion on the regime Vj 3> 1, thus dropping the corrections 
resulting from the regularization (|2.13p . This procedure 
leads to a set of constraints of the same form as in Eqs. 
(|23aj) , (|2~?j) , and ([2~gf where the terms 9 l (that in the 
classical theory are equal to CiPi) take the polymerized 
form 



sin(/i l c l ) 
®i = Pi . 
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where 
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A, according to Eq. 
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At this stage, it is worth recalling that our effective 
treatment encapsules the properties of the quantum de- 
scription from which it originates. It is thus based on 
the prescription of Refs. [23, HE] used to determine \ii. 
One can however apply directly the techniques used here 
to the corrected prescription put forward in Ref. [i~3| . 
resulting in an effective model whose relevant qualitative 
features (for our study) are in fact similar. 

Let us now rewrite the Hamiltonian constraint in terms 
of the canonical pair 
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where {bi, vj} = Sit^GSij and Vi is related with the affine 
parameter v^ defined in Eq. (12.111) via a rescaling by a 
constant. It is then straightforward to see that Qs and 
CT are the only terms in the Hamiltonian constraint that 
involve (bs, vs) and (6a-, tv), respectively; thus they are 
constants of motion. For the rest of variables, we obtain 
the effective EOM just by taking their Poisson brackets 
with the Hamiltonian constraint: 
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(3.3a) 
(3.3b) 

(3.3c) 
(3.3d) 
(3.3e) 



The evolution equations for the conjugate pair (v^ba) 
are equivalent to those of (v$,b$); one only has to make 
the replacements (<5 — > a, a — > S). The equation for is 
given by the complex conjugate of Eq. (|3.3cp . In order 
to recover the classical (noneffective) EOM, it suffices to 
realize that the dynamics dictated by the original classi- 
cal constraint (|2.8j) is given by Eqs. (|3.3j) after replacing 
sin bi with bi and cos bi with 1 (this can easily be checked 
by inspection). Remarkably, the form of the EOM for 
the inhomogeneities a m is not changed by the effective 



dynamics for each given value of the pair of constants of 
motion O5 and CT , i.e., their form is exactly the same as 
in the genuine classical case. 

At this point, it is helpful to inspect the coupling be- 
tween the different equations. The EOM for the inhomo- 
geneities and the variables {vg,be} form a system that 
is decoupled from the rest, since the variables in the ho- 
mogeneous directions 5 and a only enter these equations 
through Qg and CT . On the other hand, the equations 
for the variables corresponding to the coordinates 5 and 
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a are decoupled from each other, but coupled to the first 
system of equations through , H^ nt , and vg. In Sec. lVIII 
we will make use of this fact to solve the system defined 
by {vg, bg, a m } by their own, and analyze in this way the 
behavior of the inhomogeneities in a passage through the 
bounce. 

The system also admits another set of constants of 
motion [22j, namely, 



K„ 



(% 1 ~(X~ ri Cl_ 



(3.4) 



Quantum mechanically, each of them represents the dif- 
ference between the number of particles in the modes m 
and — m. The conservation of these quantities follows 
from the decoupling of the modes with different value of 
|m| in the dynamics and the conservation of the total 
field momentum, as required by the constraint 



Cg = ^ mK ri 



(3.5) 



m— 1 



B. Asymptotic regime 

One of the purposes of our analysis is to compare the 
properties of the universe in the distant future and past, 
regions where its evolution is expected (or at least hoped) 
to follow the standard classical trajectory. Therefore, 
with the aim at verifying this expectation, as well as to 
build suitable methods for the comparison at such distant 
epochs, it is useful to study the "low-curvature" limit of 
the EOM. This limit corresponds to small values of bg 
(for which the standard classical EOM are reproduced) 
and large vg (providing large volumes). In this limit, it 
is straightforward to solve Eq. (|3.3ap . 



vg = vg(0) exp <^ (<ds 

I 7 



(3.6) 



where vg(0) is the value of vg at a time t = 0. Substi- 
tuting this solution in Eq. (|3.3c|) . we obtain the inhomo- 
geneities: 

(3.7) 



dm — C m 6 



where c m is a complex integration constant and the phase 
Lo m depends on time trough the variable vg, 

2ttM /3AA 2 ' 3 2/3 



Introducing this solution into the definitions (|2.7[) , we 
obtain the following asymptotic behaviors, 

oo 



m— — OD 

Iml 



H Lt = T^T {' Cm ' + l C -m| COs(2w m - (f) m - 4>- m )} , 



where <p m is the phase of the corresponding complex num- 
ber c m . 



IV. BIANCHI I 

As we have already explained, Bianchi I is a partic- 
ular case of the considered Gowdy spacetimes. It can 
be described by the Hamiltonian constraint (|2.8p with 
vanishing inhomogeneities a m . Hence, the effective dy- 
namics for Bianchi I is given by the evolution equations 
p.3[) with Hq — H^ nt = 0. In this case Og is also a con- 
stant of motion, which is determined by the Hamiltonian 
constraint: 



Og = - 



(4.1) 



The system becomes thus sufficiently simple as to solve 
the EOM analytically, namely, 



= ©J cosh[a»(* - t ) + fcj 



bi (i) = 2 arctan 



(4.2a) 
(4.2b) 



with the constants fcj satisfying coshfci = l/\ sin6j(fo)|, 
rji = sgn[sin6i(i )] and 

a i = ^(6 i + e fe ), (i?j?k). (4.3) 

7 

From these solutions, it is straightforward to see that 
bounces always occur, since the hyperbolic cosine has a 
minimum. Making use of the conserved quantities, it is 
easy to obtain the value of Uj at the bounce point: 

^bounce = 2 @. ( 4 4 ) 



V. ANALYSIS OF THE BOUNCE 

In this section we focus our attention on the study of 
the bounce process. In particular, we address the follow- 
ing questions: 

(i) Does the bounce persist in the presence of inhomo- 
geneities? 

(ii) If that is the case, then is it generic? Does it hap- 
pen always or are there trajectories reaching Vi = 0, 
or at least entering the Thiemann-modified regime 
(where the quantum modifications to the inverse 
triad are important)? 

(iii) How does the presence of inhomogeneities affect the 
position of the bounce? If they happen to push the 
bounce point deeper towards the quantum regime 
(with respect to the analogous Bianchi I universe in 
vacuo), is there any lower bound on its position? 

Surprisingly, answering these questions does not require 
solving the EOM or carrying out a detailed analysis of 
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the dynamical trajectories. On the contrary, a lot of in- 
formation can be extracted just by studying the EOM, 
together with the Hamiltonian constraint, at the points 
where the necessary condition for the bounce is satisfied 
(ii = 0). These data are sufficient to define the reflec- 
tive surfaces on phase space which can never be crossed 
by any dynamical trajectory. Whenever the dynamical 
trajectory hits such a surface, a bounce or recollapse nec- 
essarily occurs. In this way, the question of whether the 
points Vi = can be reached boils down to the problem 
of elucidating whether one can find a trajectory leading 
to any of those points while avoiding the commented re- 
flective surfaces. 



A. Reflective surfaces 

Let us start with the analysis of the bounce in the 
"homogeneous" directions 6 and a. In this case, an in- 
spection of the EOM provides the answer almost straight- 
forwardly. Indeed, the evolution equation for vg (|3.3d[) 
implies that, at the reflection point (vg = 0), the function 
cos bg vanishes. Applying this to the constant of motion 
(|3.1[) and comparing with the homogeneous result (|4.4p 
(for i — 6, a) , we immediately conclude that the value of 
Vi at the bounce is exactly the same as the value for its 
homogeneous counterpart (that is, the Bianchi I solution 
with identical values of the constants of motion Qg and 

e CT ). 

This type of analysis is considerably more involved 
for the inhomogeneous direction 9. The value of vg 
at which the bounce occurs can be determined as the 
largest (positive) root of the Hamiltonian constraint, 
treated as a function of vg, upon imposing the condi- 
tion sm(be) £ {±1}- Indeed, setting this condition in the 
constraint (|2.8|) [with 6i expressed via Eq. (|3.X[> ] . and 
taking its square we get 

^(es + Qafvl = (F(v e ) - Q s Q a ) 2 , (5.1) 

where the function F(vg) is 

F(, e )^16.VG(^) 2/ %, 2/3 ^ 

, g 2 / 2 y/ 3 i4 (5 - 2) 

and contains all the dependency on the inhomogeneities. 

Given a dynamical trajectory, the variables Hq and 
Hf nt are nontrivial functions of vg. However one can as- 
sign to each trajectory the values of Hq and H^ nt at the 
bounce point. Then, on the family of all the dynami- 
cal trajectories we can introduce classes of equivalence, 
where all the trajectories with the same values of Hq 
and H^ nt at the (first, if multiple are possible) bounce 
are identified as one class. The parameter space of all 



these classes, coordinatized in particular by Hq and H^ Dt , 
is finite dimensional. On this space, the value of vg at 
the bounce will be given by the largest root of Eq. (|5.ip . 
These solutions form precisely the reflective surfaces that 
we have mentioned above. 

For the homogeneous case of Bianchi I, the function F 
vanishes and we can solve Eq. (|5 . 1 [) analytically, obtain- 
ing the result that we already know [Eq. (|4.4j) for i = 9] . 
In the most general situation, multiplying Eq. (|5.1[) by 
v^ 3 we get a fifth order polynomial P for the variable 

2/3 

v g , whose roots need to be found numerically Before 
embarking on this task, however, we note that we can 
obtain interesting qualitative results analytically. 

We first point out that there are two regions of clearly 
different qualitative behavior, depending on whether the 
function F(vg) increases or decreases the value of the 
right-hand side of Eq. (|5.1|) . Since F is positive definite, 
this depends only on the sign of the product O^O^ 

If it is negative, then the bounce (reflection) point is 
always higher than in the analogous homogeneous case. 
In particular, this result ensures that the bounce occurs 
before entering the region of small vg , where our approx- 
imation may not be valid. 

On the other hand, if Qg&cr is positive, the two terms 
on the right-hand side of Eq. (|5.1[) compete, and the 
analysis is more complicated. In order to answer the 
questions posed at the beginning of this section we have 
to analyze in detail the behavior of the roots of that 
equation. As an illustration, in Fig. [T] we plot all the 
positive roots for two particular cases (the negative roots 
are obtained from the positive ones just by a flip of sign) . 
As we will see below in more detail, the structure of the 
roots has always this form. 

Let us first coordinatize the space of equivalence classes 
of trajectories via the following reparametrization of vari- 
ables, 

tf>„ s aHl ae [0,2], Q 

where Vh — 2|0g|/3 is the position of the bounce -in 
absolute value- for the analogous Bianchi I universe in 
the 9 direction [see Eq. (|4.4jl ]. and (3 is a measure of 
the asymmetry between the two homogeneous directions. 
The values (3 = 0, 1 correspond to the degenerated cases 
when one of the homogeneous directions is suppressed 
or expanded to infinity. The bounds on a, on the other 
hand, are easily obtained from the definitions (|2.7a[) and 
(|2.7b[) (taking into account that the mode number m is 
always greater or equal than the unity). The introduced 
choice of parameters simplifies the analysis considerably, 
because the only noncompact coordinates are now Vh and 
Hi 

To arrive to conclusive results, one needs to explore the 
entire space coordinatized as above, something that can 
be done only numerically. Some analytical results can be 
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however obtained in two asymptotic regions: when Hq is 
very small and when it is very large. 

Suppose first that the inhomogeneities are small, that 
is, Hq -c 1. Then, we can neglect terms that are 
quadratic in Hq in the polynomial (15.11) . This leaves us 

2/3 

with a fifth order polynomial in v g without free term. 
In this way, we obtain one root vg = and the remaining 
(reduced) polynomial is quartic, so that we can find its 
roots via analytic methods. Apart from the zero, we ob- 
tain two complex and two real roots, out of which, given 

2/3 

the reality of v g , only the latter two arc relevant. We 
expand them as a power series in Hq. As a result, we see 
that one of them is zero up to order (Hq ) 2 and the other, 
which is the actual bounce point, is a little bit smaller 
than the position of the homogeneous bounce: 



„,bouncc 
v 9 



2/3 5/3 



v h - 16ir 2 Gr 



/3AA 



-m 



\3M J 



2/3 QsQ, 







Note that this limit of small Hq is equivalent to the case 
of large QsQa- 

In the other region, i.e., for very large Hq, it is easy to 
see that the only root of the equation tends to 



.bounce 



2M 2 / 16tt 2 7 2 G 

3 ^e 5 + e CT 



H 



(5.5) 



Returning to the general situation, the polynomial P 
can be understood as a second order polynomial in Hq. 
Finding its roots provides the values of Hq for which the 
bounce occurs at a given value of ve . These roots can be 
found analytically: 



3v, 



2/3 



Y 



-(v h ± v e ) 



(5.6) 



where we have defined 



Y:= 



\3M 



aG(6 5 + Qo 



V 2 J 



It follows that there are two positive roots for vg < vt 
and one positive root for vg > Vh . Furthermore, for large 

1/3 

vg the sole positive root is proportional to v g , as we 
have already shown via the asymptotic analysis (|5.5[) . 

Combining all these results, we get full information 
about the number of real positive roots of Eq. (|5.1[) that 
exist in different regimes. As shown in Fig.Q] one branch 
of roots starts at Vh and decreases monotonously till it 
turns around onto the branch of roots coming from zero. 
The sign of dv e /dH^ can be determined analytically in 
the following way. The derivative itself is given as the 
derivative of the implicit function P(v$, Hq) = 0, 



dve 



dP 
dHj, 



' dP 
dv e ' 



(5.7) 



and vanishes only when OP/OHq = 0. Substituting this 
condition into P(v$) = we obtain ve = as the only so- 
lution. As a consequence, in the region ve > the deriva- 
tive dvg/dHQ cannot change sign within each branch of 



roots. Therefore [and taking into account Eq. ([57 
the branch starting from Vh at Hq — is always strictly 
decreasing with Hq. 

By the very same method one can also prove that the 
line F(vg) = <dsO a [with F given in Eq. (|5.2|) ] never 
intersects any of the roots of Eq. (|5.1|) outside the point 
0,v g = 0. 



Hi 



B. Minimum bounce point 

In the previous subsection we gathered sufficient data 
to characterize the reflective surfaces on phase space. In 
particular, in the case Os®a < 0, we proved that the 
bounce point is always larger than its homogeneous coun- 
terpart. On the other hand, for 0,58,7 > 0, and iden- 
tifying in the same equivalence class all the dynamical 
trajectories with identical values of Hq and Hf nt at the 
(candidate) bounce point, we showed that the problem 
of finding a lower bound for the bounce can be recast as 
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FIG. 2. The ratio between the minimal root and the value of 
the fixed homogeneous bounce Vh — 10 5 in terms of a and /3. 



that of determining the infimum of the largest real solu- 
tion to Eq. (15.11) . Now, have this largest solution been 
a continuous function of Hq and Hf Dt , its infimum would 
have provided exactly the lower bound on the bounce po- 
sition for the dynamical trajectories. However, the func- 
tion presents a discontinuity, a fact which complicates 
the situation. 

Fortunately, an extensive numerical analysis (de- 
scribed in detail in Sec. IVI[) shows that, on any given 
dynamical trajectory, the relative change of Hq is quite 
small. As a consequence, the trajectories can be approx- 
imated by vertical lines on the plane Hq-vs- Moreover, 
the roots of Eq. (|5.1[) change slowly with H^ nt . In this sit- 
uation, a lower bound on the position of the (first) bounce 
(if multiple bounces are possible) for fixed Vh and /3 is in- 
deed given by the corresponding infimum of the largest 
solution to Eq. (|5.ip . Owing to our approximations, 
nonetheless, the equality between the former bound and 
this infimum is also approximate only. However, it is 
not necessary to determine the bound exactly. Its order 
of magnitude is enough to establish the kind of qualita- 
tive results we are interested in (furthermore, our result 
will be checked in Sec. IVII Bl by performing a numerical 
analysis of large populations of dynamical trajectories) . 

As we have seen, for small Hq the bounce point is 
smaller than its homogeneous counterpart, and this value 
decreases with Hq. This provides the largest root of Eq. 
(|5.ip until the number of positive roots changes. Be- 
yond this discontinuity, the new largest root belongs to 
the branch that starts at zero for vanishing Hq (as can 
be checked in Fig. [T]), and becomes an increasing func- 
tion of the inhomogeneities Hq. Therefore, the reflecting 
boundary exists for all values of Hq, and the infimum 
of the desired function occurs exactly at the commented 
discontinuity. We accordingly adapt our search methods 
to these considerations. 

On a general level, the discussion reduces to study the 
function vs/vh, where vb is the largest root of Eq. (|5.ip 
for fixed parameters (Hq, a, /3, Vh). The problem is split 
into several steps as follows. First we fix the values of 
all the parameters except Hq and find the minimum (or 
rather infimum, because it happens at a discontinuity) of 
vb/vh with respect to Hq, which we call v m i n . Next we 



a 
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a 

FIG. 3. The ratio Vmi n /vh for a fixed value of /3 = 0.5. In the 
first plot, Vh appears in a logarithmic scale, but in a linear 
scale in the others. The second and third plots are in fact the 
same, only the perspective is modified. 

analyze the properties of v m i n (a, f3,Vh). Since this anal- 
ysis must be performed numerically, we choose natural 
units with G = 1 and l p = 1 , whereas for the Immirzi pa- 
rameter we take the value deduced from the computation 
of black hole entropy H3, i.e., 7 = 0.23753295 . . . 

In order to find the value of Hq at which w m ; n is at- 
tained, we first note that it must be smaller than Hq, 
which is defined as the maximum value of Hq allowed by 
the implicit function 

F(v e ,Hl) = Q s e a (5.8) 

for positive vg. This can be clearly seen in Fig. Q] (see 
also the discussion at the end of Sec. IV Al for more pre- 
cise conclusions). Therefore, to find f m in it is enough to 
explore numerically the interval Hq G [0, Hq]. 

This numerical analysis was performed in the following 
way. The domain H^ G [0, H$] was divided in 1000 uni- 



11 




O.ll 

0.05* 



5 10 15 20 25 hl v 

FIG. 5. The ratio Vmin/vh for fixed values of a = 0.1 and 
f3 = 0.5. This is a section of the previous plots which shows 
the infimum of « m i n ~ 0.05v h . 

form intervals and the largest root of Eq. (I5.1|) was found 
at each point ni^/lOOO, for n = 0, . . . , 1000. The min- 
imum of those roots was chosen, with its corresponding 
value of Hq = i?™ ln , and the study was then repeated 
in the interval [H^ in - i^/lOOO, Hj? in + i?o/1000]. In 
Figs. [U 01 131 and [5] we show the behavior of the function 
v m in{&, ft, Vh)/vh with respect to the different parame- 
ters. In each plot, one of the three parameters (a,f3,Vh) 
is fixed. 

The main result found is that the value of the minimum 
bounce point is bounded from below by t> m in « 0.05u/,., 
as it is evident from Fig. [5] This bound is approached 
at values of a and Vh that are small. In fact for large 
Vh and/or a, the ratio v m i n /vh tends to a number close 



to 0.15. In conclusion, we have shown for the inhomoge- 
neous case that vg at the bounce is always larger than a 
5% of its homogeneous counterpart, or in other words, it 
is never two or more orders of magnitude smaller. 

Let us emphasize that the analysis performed in this 
section does not actually solve the EOM. Instead, we only 
provided reflective boundaries on phase space. Nonethe- 
less, since these boundaries are discontinuous, we can not 
exclude that a dynamical trajectory may slip between dif- 
ferent branches of roots and reach the singularity vg = 
in that way. However, this possibility would require an 
extreme fine tuning. Furthermore, it is worth noting that 
such a possibility would require that the value of Hq on 
the trajectory decay to zero in the evolution. This implies 
that, in any case, the system would need to homogenize 
itself in order to reach vg = 0. This behavior might then 
be interpreted as a particular realization of the BKL con- 
jecture [38] , restricted to the case of the Gowdy spacetime 
investigated in this article. 

Moreover, the above situation is possible only on tra- 
jectories for which all constants of motion K m [defined in 
Eq. (|3.4[) ] vanish. An example of a near-critical dynami- 
cal trajectory of the system can be found in Fig. [7j where 
the initial data are chosen so that the bounce happens 
near the discontinuity in the largest root of Eq. (|5.ip . 
In this case, the system bounces and recollapses back 
and forth several times before escaping to the region of 
asymptotically large vg. 



VI. NUMERICAL STUDY OF THE DYNAMICS: 
THE BOUNCE 

Despite the fact that a lot of information about the 
system can be extracted just from the Hamiltonian con- 
straint and the EOM without integrating them, a full 
answer to most of the questions posed at the beginning 
of this paper requires the numerical determination of the 
dynamical trajectories. In this section we present the nu- 
merical techniques that we have used to solve the EOM, 
and display and discuss explicit examples of solutions. 
Here, we focus our attention on the behavior of the ho- 
mogeneous variables (namely Vi). The analysis of the be- 
havior of the inhomogeneity modes is postponed to the 
next section, since it requires more specific methods. 

Let us start with the initial value problem formulation 
and the description of the integration method. 



A. The initial value formulation 

The analysis of the asymptotic behavior of the solu- 
tions of Sec. IIIIBI revealed that, in the large volume 
regime, each of the inhomogeneities a m becomes an os- 
cillating function of constant amplitude, with a time- 
dependent frequency that is proportional to v 2 J 3 (13.71) . 
Since one enters that regime when the solutions to Eq. 
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(|3.3|) are constructed, these oscillations require elabo- 
rated integration methods, forcing one to refine the inte- 
gration step and increasing the computational time sig- 
nificantly. To avoid this difficulty, we implement two 
changes in the formulation of the system. First, we 
change the time t to a new time t' defined as, 

d _ 1 d 

dt' \v e \ 2 / 3 dt' [ ' ' 

In practical terms, this is easily done by dividing the 
right-hand side of Eqs. (|3.3p by |f;e| 2 ^ 3 - With the new 
choice of time the asymptotic behavior of the inhomo- 
geneity modes takes the following form: 

a m = c m e-™™ t \ (6.2) 

where the frequency u) m is 

(3M\ 2 / 3 

V) • (6 - 3) 

This frequency depends on the absolute value of m but 
not on time. As a result, now each inhomogeneity mode 
oscillates with a constant frequency. 

Next, it is convenient to introduce a change of vari- 
ables from a m (t') to c m (t') := a m (t')e tulmt . The asymp- 
totic relation (|6.2I) implies then that the new variables 
c m (t') approach asymptotically the constants c m . This 
allows one to significantly improve the efficiency of the 
numerical integration, because the oscillatory behavior is 
removed in the regime of large v, permitting one to adopt 
larger integration steps there. 

As we explained in Sec. II, the system of equations 
that we want to solve numerically has several conserved 
quantities: the Hamiltonian and momentum constraints 
[ (|2.8I) and Q3.5P respectively], the constants K m [given 
in Eq. (|3.4[) ]. Og, and Oo-. In principle, one can solve 
the constraints they impose on the system, decreasing 
in this way the number of variables. Here, however, we 
use a free-evolution scheme and employ the conservation 
of these constants of motion as a (quantitative) test on 
the correctness of our numerical method. Therefore, the 
equations that we want to integrate numerically are those 
for the variables {vi,bi} [ (|3.3ati3.3b|) and (|3.3dH3.3c"l) ] in 
terms of the new time variable t' . In addition, the equa- 
tions (I3.3cp are replaced with the evolution equations of 
the new inhomogeneous variables c m , which are 

dc m iG f 2 \ 2/3 

dt' 7 2 HK| 4 /3 \3MJ ( 6 .4) 

x (e 5 + e CT ) 2 (c m + c*_ m e 2iat '). 

To compare the solutions of these effective equations 
with their counterpart in General Relativity, we also need 
to solve the classical EOM of our Gowdy spacetimes. In 
practice, these can be obtained from the effective equa- 
tions introduced above just by making the transforma- 
tions {sin&i — > bi, cos bi — ¥ 1} (see the discussion in 
Sec. IIII Al for the validity of this method). 



As initial conditions for our simulations, we select data 
corresponding to classical universes. In particular, we al- 
ways start the simulation at large values of the variables 
Vi and small values of the variables bi . The crucial factor 
responsible for the variety of studied solutions is found in 
the initial conditions for the inhomogeneity modes. On 
the one hand, different cases come from a different num- 
ber of excited modes. On the other hand, the amplitude 
of the modes can also vary from case to case. 

The constraints (|2.8[) and (|3.5p must be satisfied, in 
particular, on the initial slice; hence they must be taken 
into account in the construction of admissible sets of ini- 
tial data. The exact construction is as follows. First, 
we build data with just n excited modes, that is, the 
only nonvanishing modes c m are those with \m\ < n, for 
some fixed number n. Except for Im(c„), the values of 
all the initial data c m are constructed via a random num- 
ber generator, with a probability density function given 
by a Gaussian of vanishing mean value. Finally, we ob- 
tain the remaining piece of data, Im(c„), by solving the 
momentum constraint (I3.5p . 

The amplitudes of these modes are of the order of the 
standard deviation of the probability distribution func- 
tion, which is the remaining freedom in the construction. 
We exploit this freedom, changing this deviation in sev- 
eral ways (see Appendix [A] for details) in order to check 
the robustness of the result. With the same aim, in some 
simulations we actually use the flat distribution instead 
of the Gaussian one, but this choice does not lead to any 
relevant difference in the results. 

Given the initial data for the inhomogeneities obtained 
in this manner, we next fix the remaining free parame- 
ters: (v$,b$), {v a ,b a ), and vg. The initial value of the 
remaining variable bg is then determined by the Hamil- 
tonian constraint (|2.8I) . 

It is important to recall that we arc interested only 
in initial data which provide universes in the classical 
evolution regime. Therefore, after constructing the data, 
we have to verify that Vj&j ~ Vi sin 6,; for all directions i, 
up to the desired accuracy (this is a condition additional 
to the ones discussed above). 

Since the universe bounces in the effective trajectory, 
in order to compare it with its classical counterpart we 
have to solve the classical EOM for two cases. Namely, 
the case with the same initial data and the case where 
those data are provided by the effective trajectory at 
large distant times in the future. This latter case gives 
a good approximation to the classical universe that the 
effective trajectory approaches after the bounce. These 
two classical universes are related by a flip of sign in the 
quantities 85 and Q a . 

Once the initial data are selected, the system of EOM 
is integrated via the built-in Mathematica adaptive func- 
tion NDSolve, which by default uses an Adams method 
and a Gear backward difference method, switching dy- 
namically between them [39]. For our simulations, specif- 
ically, we set the error control options to AccuracyGoal— >■ 
10, PrecisionGoal— > 12, and WorkingPrecision— > 25. 



13 



Since only a finite number of modes has been excited 
initially and the linear equations (16.4)) are homogeneous 
within each set of modes with identical values of |m|, only 
a finite number of EOM needs to be integrated. 

Let us end this subsection with a remark. In order to 
overcome technical difficulties (infinite number of EOM) 
only a finite number of modes has been allowed to be 
excited. One could then argue that the system has been 
effectively truncated to one with only a finite number of 
degrees of freedom, and thus its behavior might be signif- 
icantly different from its field-like counterpart. Even so, 
and at least from our effective perspective, both systems 
(the truncated one and the infinite dimensional system) 
should behave similarly. This is supported by the fol- 
lowing argument. Since the value of Hq provided on the 
initial section is finite, the series (|2.7al) which defines it 
must converge at least near the initial time. This im- 
plies that there exists a finite positive integer I such that 
\a m \ > |»m+x| f° r an l TO l ^ On the other hand, the 
form of the EOM ensures that the growth of each individ- 
ual mode cannot be faster than exponential. Hence the 
quantities Hq and may diverge only if the amount 
of the relative change of | a m | across the bounce increases 
sufficiently fast with \m\. However, as we will show in 
Sec. I VII B[ statistically there is no detectable evidence 
that this change depends on the value of \m\ at all. 

Furthermore, the inhomogeneities interact only within 
pairs (m, — m), and their amplitude can be "pumped up" 
only through the backreaction of the homogeneous back- 
ground. As we will see in Sec. IVII1 this process occurs 
only near the bounce and does not depend on the number 
of excited modes (as the EOM for the homogeneous vari- 
ables involve only the total energies Hq and Hf nt ). As a 
consequence, one can extrapolate the commented results 
about the dependence of the amplification on |m| to the 
case of an infinite number of modes. This indicates that 
a change in the mode energy distribution that is suffi- 
ciently strong as to destroy the convergence of Hq occurs 
with vanishing probability. Therefore, the contribution 
to the system dynamics of all the modes with \m\ larger 
than a certain threshold ought to be negligible. 



B. The dynamical bounce 

As a way to present the results of the numerical anal- 
ysis of the dynamical trajectories, let us give an example 
of a typical solution and discuss its behavior. 

In Fig. |5] we show a generic example of regular evo- 
lution in a logarithmic plot on the (vg,vs) and (vg,v a ) 
planes of phase space. This specific plot represents a 
system with five modes. The Gaussian width a of the 
random number generator that provides us with the ini- 
tial conditions for the inhomogeneities has been cho- 
sen equal to 10 _1 , and leads to the approximated val- 
ues H$ « 0.3846ft and H mt « 0.02719ft (this makes 
a sa 0.0707). We observe that the effective trajectory, 



represented by a black (continuous) line, asymptotically 
connects two different classical trajectories. Namely, 
the effective trajectory converges to two distinct classi- 
cal ones in its distant past and future, respectively. As 
we have already mentioned, these two classical trajecto- 
ries are related by a flip of sign in Q$ and CT which, 
for this particular example, take the values &s = 300lp 
and CT = 100/p. The evolution begins at vg = I0 6 lp, 
Vs = 10 4 ^, and v a = 10 4 /p. In the (vg,vs) plane, the ef- 
fective trajectory follows a classical one corresponding to 
a contracting universe almost until a bounce in vg hap- 
pens. Near this bounce, the solution departs from the 
classical (diverging) trajectory and, very quickly after vs 
bounces, approaches the new expanding classical trajec- 
tory. 

The numerical simulations have been performed for 
various values of the inhomogeneities, with data gener- 
ated by a Gaussian random number generator of chang- 
ing width. This leads to a population of different tra- 
jectories in phase space, for which the behavior of the 
bounce agrees with the results of the analysis of Sec. fVl 
With exception of the near-critical trajectories, the be- 
havior is qualitatively the same as described for the ex- 
ample presented in Fig. [5] 

To complete the analysis of this section, let us consider 
now a near-critical solution. An example of such a kind of 
solution is shown in Fig. [7] In that particular simulation, 
we have chosen initial data with K m — 0, a necessary 
condition for the system to be potentially able to reach 
the singularity. The plot shows a zoom of the dynam- 
ical trajectory, presenting the variation of Hq in terms 
of vg near the mentioned point. During the evolution, 
the system becomes trapped between the two "leaves" 
of the reflective surface. It then repeatedly bounces and 
recollapses, upon hitting the lower and upper leaves cor- 
responding to the different branches of roots of Eq. (|5.ip . 
After several of these cycles, the system is finally able to 
escape from this region to the sector of classical regime. 
It is extremely difficult to say something general and pre- 
cise about the evolution of the system near this "throat" , 
formed by the distinct branches of roots. Actually, this 
evolution is highly sensitive even to very small changes of 
the mean value of the inhomogeneities during the bounce. 
In principle, the possibility exists that the trajectory may 
even go all the way "down the throat" and reach the 
singularity in an infinite sequence of bounces and recol- 
lapses. However, such critical trajectories require a fine 
tuning of parameters; thus (apart from being unphysi- 
cal due to the own limitations of the effective treatment) 
they should be disregarded based on statistical consider- 
ations. 



VII. NUMERICAL STUDY OF THE 
DYNAMICS: THE INHOMOGENEITIES 

Once we have established that the bounce persists, it is 
natural to ask how the structure of the inhomogeneities 
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FIG. 6. The black (continuous) curves are dynamical effective trajectories which, after bouncing in all the three directions, 
converge to the classical trajectories displayed in (dashed) red. The initial data are vg = 10 6 ll, v s = 10 4 ^, and v a = 10 4 ^, 
with constants of motion = 300ip and B CT = lOOip. Only five inhomogeneity modes have been excited. In order to obtain 
initial conditions for those inhomogeneities, the width of the Gaussian random number generator has been chosen to be 1CP 1 . 
This has led to the initial values Hq « 0.3846ft and H lnt ~ 0.02719/i (a « 0.0707). The two classical trajectories are related by 
a flip of sign in the constants <ds and O ct . 
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FIG. 7. An example of the near-critical behavior of the evo- 
lution near the throat formed by different branches of roots. 
The system bounces and recollapses several times before going 
back to the region with large values of vg. 

changes across it. This section explores this issue. In 
particular, we will focus our discussion on the amplitudes 
\a m \ of the different inhomogeneity modes, since these 
quantities contain physical information about the energy 
of those gravitational waves modes. 

A. General description 

Let us first study the behavior of the amplitude of a sin- 
gle inhomogeneity mode. A typical dynamical trajectory 
is displayed in Fig. [FJ the evolution being represented 
with respect to v$. As one can see, before and after 
the bounce (in ve) the trajectory of \a m \ describes si- 
nusoidal oscillations about a constant (conserved) mean 
value. The amplitude of these oscillations behaves like 
Vg e for certain positive number e. In particular, for 
large values of vg, \a m \ tends to a constant value |c m |, 
as proves the asymptotic solution (|3.7j) . On the other 
hand, the amplitude of the oscillations is amplified as 
the trajectory approaches the bounce, reaching the max- 
imum there. At that point, the oscillatory behavior is 
interrupted temporarily. This may result in a change 
of the mean value of the considered mode. In particu- 



FIG. 8. Typical evolution of the amplitude \a m \ of two modes 
of the inhomogeneities. During the evolution, away from the 
bounce, its mean value (obtained by subtracting the oscilla- 
tions) is conserved. However, we see two distinct types of 
behavior for the change of this mean value trough the bounce 
in vg. In case (a), the mean value changes by an amount 
which is almost one amplitude of the oscillations, whereas it 
remains practically unchanged in case (b). 

lar, Fig. [5] shows the two different possible behaviors of 
the modes when crossing the bounce. In this example, 
the mean value increases in one case by a full amplitude 
of the oscillations, whereas in the other case the mode 
traces back its pre-bounce trajectory almost perfectly af- 
ter the bounce, and thus its mean value is approximately 
preserved. 

In order to understand the physical consequences of 
these different behaviors, we have integrated the evolu- 
tion for a set of initial data where only the phase of a 
given mode is allowed to vary, namely, by a phase shift 
tp. We have then measured the change of the mean value 
of this mode after the bounce, treated as a function of 
the phase shift. Technically, the mean value has been 
identified with the value at the point where the second 
derivative of the mode amplitude \a m \ vanishes. 

This study again reveals two distinct types of behavior, 
presented in Fig. [§J In the first case, that corresponds 
to the situation where the oscillations of \a m \ near the 
bounce are small in comparison to its mean value, the 
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change in this mean value is antisymmetric under a ro- 
tation of the initial phase shift, ip — > tp + tt. Hence, for 
large samples, the mean value of the amplitude of the in- 
homogeneities, and thus of Hq, is conserved through the 
bounce. The second case corresponds to the opposite sit- 
uation. The sinusoidal function is antisymmetric around 
certain positive axis; thus it has a positive integral. This 
implies that, in this case, the inhomogeneities are statis- 
tically amplified. We will further investigate this effect 
in the next section, carrying out a systematic statistical 
analysis. 



B. Statistical analysis of the energy distribution 

In modern cosmology, a great deal of attention is paid 
to the analysis of the energy distribution of the inho- 
mogeneities, with the phases of the modes being often 
ignored because, essentially, they do not affect the obser- 
vations. Furthermore, in the class of cosmological mod- 
els considered here, those phases influence the universe 
dynamics only through the quantity if-L, which is sup- 
pressed with respect to Hq (which does not depend on 

the phases) by a factor v^J 3 . As a consequence, in the 
far future and past, where the universe is in a classical 
regime, their effect is negligible. 

One has to recall that they are still part of the initial 
data (see Sec. IVI Ap . so they are necessary to determine 
the dynamical trajectory in a unique way. However, it 
is natural to ask what can be said about the evolution 
once they are ignored. In particular, given data in the 
distant past determining the energies of the gravitational 
wave modes (encoded just in |a m |), one can ask about 
their value and properties in the distant future. Since we 
are neglecting part of the initial data, the nature of the 
analysis and of the results necessarily becomes statistical. 

With this motivation in mind, in this section we are 
going to analyze the statistically averaged change of en- 
ergy of the individual modes during the transition from 
the distant past to the distant future. This is a basic 
piece of information, necessary to assess the changes in 
the energy distributions. 

To achieve this goal, we employ a Monte-Carlo anal- 
ysis, that is, we carry out a large number of numerical 
integrations of the system using different random initial 
data in order to cover a portion of phase space as large 
as possible. Since the statistical behavior has proven not 
to depend on the number of excited modes (as one can 
check by inspection), we have performed the bulk of the 
simulations just for the case in which only the first five 
modes are excited, i.e., we assume a m = for all |m| > 5. 
The results of these simulations, nonetheless, have been 
verified on much smaller samples of data but with already 
twenty modes excited. 

In addition, to ensure the robustness of the results, the 
simulations have been performed for six different forms 
of the energy distributions. For each individual trajec- 



tory, the construction of the set of initial data has been 
specified along the lines explained in Sec. IVI Al In more 
detail: 

• The value of the constants of motion 0$ and CT have 
been given in terms of the partly compactified coordi- 
nates Vh and (3, as defined in Eq. (|5.3I) . and whose values 
have been chosen by a random number generator with flat 
distribution. In particular, the value of the anisotropy 
factor /3 has been selected from the interval [0.1,0.9] to 
exclude the uninteresting near-degenerate cases from the 
analysis. The value of Vh has been chosen as the expo- 
nential of a number (again provided by a random number 
generator of flat distribution) from certain interval. This 
interval depends on the data set (see the specification of 
the inhomogeneities below) and is equal to: [0,logl0 9 ] 
for cases (i)-(iv), and [0,logl0 5 ] for (v)-(vi). 

• The initial values of the inhomogeneities [all but 
Im(a5), which is determined by the diffcomorphism con- 
straint] have been obtained (real and imaginary part sep- 
arately) using a random number generator with a Gaus- 
sian distribution of vanishing mean value. The width a 
of this Gaussian has been selected with the use of six dif- 
ferent algorithms, which in fact define our data sets. For 
each of these algorithms, a has been specified as follows: 
(i) a is a random number (again with flat distribution) 

from the interval [0, 1000], (ii) a = 100, (Hi) a = 0.1^ /3 , 

(iv) a = 0.5vf/ 3 , (v) a — 2vf/ 3 , and (vi) a = Vh- 

• Given the above data, we have solved Eq. (|5.ip to 
estimate the bounce point, choosing the initial value of v g 
as 10 7 times the obtained value. In this way, it is assured 
that the value of vg is large enough as to guarantee that 
the initial point lays in the classical regime. 

In the last four of the above sets of initial data for 
the inhomogeneities, we have related the value of a with 
the homogeneous bounce point vt- The reason for this 
is that, in a universe with large Vh, the inhomogeneities 
should be larger in order to have a similar effect in its dy- 
namics. This is motivated by the fact that the amplitude 
of the inhomogeneities (or equivalently of Hq) by itself is 
not a good measure of the departure of a universe from 
homogeneity. For that purpose one should rather con- 
sider the relative values of the terms that contain the ho- 
mogeneous and the inhomogeneous contributions in the 
Hamiltonian constraint (|2.8p . 

The results corresponding to the six different initial 
data selections are compared in Appendix [Al 

Once we have specified all the initial data, the dynam- 
ical trajectories have been integrated until a final time 
tfin defined by the relation vg(tn n ) — vg(to). Since our 
primary interest is the behavior of the inhomogeneities, 
only the subset of EOM for {vg, bg, a m } has been solved. 
This is possible because of the decoupling explained in 
Subsec. IhTAI 

We have integrated 400 dynamical trajectories in case 
(i) and 1000 trajectories in each of the cases (ii) — (vi), 
thus reaching the total number of 5400 trajectories. In 
order to test the correctness of the numerical evaluation, 
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FIG. 9. The dependence of the change of the amplitude through the bounce on the initial phase shift of the mode. 



we have checked the conservation of the Hamiltonian and 
momentum constraints, as well as of the constants K m 
(this check is done by comparing their values at the initial 
and final times, to and ifi n ). Imposing this requirement 
with an absolute precision of 10 turns out to rule out 
831 trajectories, that we consider incorrect, thus leaving 
us with a final population of 4569 dynamical trajectories 
to analyze. 

For each element of this population, we have computed 
the relative change, 

a I | |fflm|fan) ~ \a m \(t ) 

\a m \[to) 

of the inhomogcneity amplitude \a m \, analyzed sepa- 
rately for each excited mode. The results have then been 
compared with the ratio between the actual bounce point 
Vb [not to be confused with the solution of Eq. (|5.1[) . vb] 
and its homogeneous counterpart. This ratio, following 
from the analysis of Sec. lVB] provides us with a measure 
of how inhomogeneous a given universe is. 

The final results are displayed in Figs. [TUl and fTTl where 
A | a m | is plotted against Vb/vh for the whole population 
of trajectories and in each of the cases (i)-(vi), corre- 
sponding to the different algorithms used to generate de 
initial data for the inhomogeneities. Looking at Fig. [TOl 
one can distinguish two distinct regimes: 

(a) inhomogeneity dominated (large Vb/vh), for which the 
dynamics around the bounce is dominated by the 
content of gravitational waves, and 

(b) near-homogeneous (small Vb/Vh), for which those 
waves introduce only small corrections to the vacuum 
Bianchi I dynamics. 

These two regions are clearly separated by an intermedi- 
ate region where A|a m | drops sharply. 

In the first of these regimes, the relative changes are 
very small; thus, by the results of Sec. IVII Al the data 
points are in the regime where the average of A|a m | over 
the dependence on initial phases is ensured to vanish. On 
the other hand, in the second regime above, we see that 
a large fraction of the data points exceed A|et m | = 1, and 



hence they correspond to the case where the antisymmet- 
ric behavior of Al function of the phase differ- 
ences is generically lost, and its average becomes strictly 
positive. Therefore, in the near-homogeneous regime the 
quantum geometry effects around the bounce enhance the 
inhomogeneities. 

An intuitive explanation of this mechanism comes from 
the analysis of the previous subsection. As far as we have 
been able to observe, the maximum change that a given 
mode \a m \ can experience at the bounce is of the or- 
der of the amplitude of its superimposed "oscillations" 
in the nearby region. This latter amplitude increases 
like Vg e , and nearly homogeneous universes have much 
smaller value of vg at the bounce than the highly in- 
homogeneous spacetimes. This allows for larger relative 
changes of \a m \. When the amplitude of the oscillations 
in Fig. [5] exceeds the mean value of \a m \, the antisym- 
metry in the dependence on the phase difference is lost, 
because \a m \ must always remain positive by definition. 
The breaking of this antisymmetry opens the possibility 
of an amplification. 

Let us finally remark that the minimum bounce point 
that we have observed in the studied population of dy- 
namical trajectories is Vb — 0.070u/j- This value is con- 
sistent with the lower bound given in Sec. IV Bl for the 
bounce point and, in fact, it is quite close to it. 



VIII. SUMMARY AND DISCUSSION 

We have studied a class of cosmological inhomoge- 
neous spacetimes -the vacuum Gowdy universes of 3- 
torus topology- taking as starting point a first order ef- 
fective theory, built out of a quantum model that has 
been obtained by implementing the hybrid quantization 
scheme of LQC. It is worth emphasizing that the pre- 
sented treatment constitutes the first nonperturbative ef- 
fective description of an inhomogeneous spacetime within 
the loop framework which has been constructed in a sys- 
tematic and controllable way out of a rigorously defined 
quantum system. Using this description, we have inves- 
tigated the dynamics of a large population of universes, 
focusing on two issues. 
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FIG. 10. In these plots we show, in a linear and a logarithmic scale respectively, the absolute and the relative change in 
the amplitude of the inhomogeneities \a m \ produced when the bounce is crossed, represented against the ratio between the 
actual bounce point «(, and its homogeneous counterpart Vh,. We note that, for small Vb/vh (almost homogeneous universe), 
A|o m | can exceed the unity, regime which corresponds to the situation when the inhomogeneities are amplified. On the other 
hand, for large Vb/vh (very inhomogeneous universe), |A|a m || < 1CP 4 , and thus the amplitude of the inhomogeneity mode is 
approximately conserved. The horizontal line in the right plot corresponds to the data points where |A|a m || < 10 -16 , and 
could not be measured accurately owing to floating point precision limits. 



The first of them is the verification of the persistence of 
the bounce phenomenon observed in homogeneous LQC 
in spite of the presence of inhomogeneities. With that 
aim, we have compared the properties of an inhomoge- 
neous universe with those of its Bianchi I counterpart, 
namely, the homogeneous background on which the grav- 
itational waves responsible for the inhomogeneities prop- 
agate. The analysis of the EOM and the constraints 
have proven that the bounce is indeed present, perhaps 
with the possible exception of certain critical trajecto- 
ries. In two of the directions in configuration space (de- 
fined in terms of the areas of large 2-tori) the bounce 
of the Gowdy universe and its Bianchi I counterpart 
coincide. For the remaining ("inhomogeneous") direc- 
tion, the exact results depend on the considered region 
of phase space. In particular, in the region where the 
constants of motion Qg and 6 CT [see Eq. (|3.1[) ] satisfy the 
inequality <d$O a < 0, the bounce always happens further 
away from the analog of the classical singularity than in 
the Bianchi I counterpart. In the other case, OgQo > 0, 
the bounce can actually happen closer to the singularity, 
and determining its properties has required a combined 
analytical and numerical analysis. We have shown then 
that there exists a lower bound on the ratio between the 
bounce positions in our system and in its Bianchi I coun- 
terpart. In terms of the parameters vg (|3. 21) . this bound 
is approximately 0.05, and is valid in the whole space of 
solutions, except possibly for an extremely small region 
around the zero measure set of critical trajectories. 

Let us now comment on those critical solutions. They 
represent a very special and quite peculiar class of tra- 
jectories, as the universe following them actually can col- 
lapse all the way down to a singularity through a very 
intricate, infinite sequence of bounces and recollapses. 
However, owing to the intrinsic limitations of the effec- 
tive treatment, they cannot be considered as physical. 
Indeed, they pass through the bifurcation points, where 



minute changes of initial data can modify the global tra- 
jectory drastically. In consequence, the effects of the 
higher-order state dependent parameters, here neglected 
in the effective description, can play a significant role in 
theses cases, and the effective description itself ceases to 
be valid. 

Once the presence of the bounce as a generic feature 
has been positively confirmed in the considered model, we 
have addressed the second issue, namely, whether there 
exist significant changes to the structure of the inhomo- 
geneities across the bounce. Here, our investigation has 
been mainly focused on comparing the distributions of 
the amplitudes of the inhomogeneity modes (which en- 
code the information about the distribution of the energy 
of the gravitational waves) in the distant future and past 
of the universe. To discuss this question, we have con- 
verted the deterministic dynamical system into a statis- 
tical one by averaging over phases in each of the modes. 
In order to study the properties of this statistical system, 
the Monte-Carlo method has been applied. First, a large 
population of dynamical trajectories has been calculated 
numerically. And then, the data representing the asymp- 
totic behavior of the inhomogeneity modes has been ex- 
tracted and used to probe a statistical relation between 
the relative change of the amplitude of the modes and 
the measure of the departure from homogeneity. 

Our analysis has revealed the existence of two separate 
regions in the space of solutions. One of them, denoted 
as inhomogeneity-dominated, contains those Gowdy uni- 
verses which bounce at volumes much greater than their 
homogeneous counterparts, and for which the bounce 
process itself is driven by the inhomogeneities. In that 
regime, the amplitude of the modes is statistically pre- 
served. On the other hand, in the near-homogeneous re- 
gion, where the inhomogeneities produce only small cor- 
rections to the evolution of the homogeneous background, 
the inhomogeneities turn out to be statistically amplified. 



18 



The mechanism of this amplification has been understood 
by analyzing the behavior of individual modes. Remark- 
ably, in our statistical investigation, we have not detected 
any measurable difference that could be assigned to a de- 
pendence on the mode number of the inhomogeneity. 

The above results provide a potentially interesting 
mechanism for structure generation. If the universe is 
sufficiently homogeneous, the bounce process enforces the 
growth of inhomogeneities, in this case pumping energy 
to the gravitational waves. Only when those inhomo- 
geneities reach a sufficiently large amplitude, the growth 
stops and their structure stabilizes. 

The very same mechanism allows one also to shed 
some light on the process of entropy growth through the 
bounce. Indeed, in the system under study, the energy of 
the inhomogeneities provides an intuitive entropy mea- 
sure. The results of the above paragraph can be thus 
reinterpreted in the following way. In the process of the 
bounce, the entropy of the universe would increase till 
it reaches certain critical order of magnitude. Then, the 
bounce process would preserve it. Of course, this sta- 
tistical mechanism is applicable only for populations of 
universes, not for a single one. However, it may become 
feasible in the case of classically recollapsing models, like 
for example Gowdy universes of the topology S* 1 x S 2 
and 1S 3 [40] - In these cases, the infinite chain of bounces 
and recollapses (within the evolution of one universe) can 
make the statistical process physically meaningful. 

All these results, though very promising, should be 
treated only as preliminary for various reasons, that we 
now comment. 

First, our effective treatment is only a first-order the- 
ory, where all the state-dependent quantities have sim- 
ply been dropped. Therefore, its applicability and re- 
liability needs to be further confirmed, either by a sys- 
tematic derivation of the effective dynamics with careful 
estimates of higher-order corrections or by a comparison 
with the genuine quantum evolution in a domain suffi- 
ciently large as to be representative. We already know 
from our discussion in this section that the treatment is 
to fail in the near-critical regime. The verification be- 
comes particularly important, since we are considering 
a system with an infinite number of degrees of freedom, 
and it is only in a series of systems with a finite num- 
ber of them that it has been possible to verify so far the 
correctness of the effective dynamics. Fortunately, in the 
considered model, the structure of the dynamics of the 
inhomogeneous degrees of freedom is sufficiently simple 
(modes couple only in pairs) to believe that it is possible 
to extrapolate the results about the validity of the effec- 
tive description from the quantum mechanical systems. 
Furthermore, the decoupling between (pairs of) modes 
opens the possibility of verifying this validity by apply- 
ing technically manageable inductive methods, based for 
example on a combined effective/genuine treatment [4l|. 
or adapting the systematic effective treatment of Ref. 
This will be the subject of future investigations. 

The above limitations are connected with an impor- 



tant and potentially relevant physical problem: the zero 
point energy. Namely, as discussed in Sec. Ill BL the in- 
homogeneity modes form a Fock space and each of them 
(after having implemented a convenient conformal trans- 
formation) behaves like a scalar field with a time depen- 
dent mass, with a ground state that is not dynamically 
invariant. As a consequence, one may expect the (dy- 
namical) generation of a nonzero energy density, which 
could radically affect the system once all the modes are 
taken into account. This effect, certainly happening in 
some geometrodynamical systems when one adopts cer- 
tain choices for the Fock quantization of the inhomo- 
geneities, is currently under investigation. Concerning 
the choice of Fock quantization, it is worth emphasizing 
again that the quantization adopted here is the only one 
that ensures a unitary quantum evolution -without (at 
least some) divergences- after a suitable deparametriza- 
tion, which results from a convenient choice of internal 
time. On the other hand, introducing a polymer quanti- 
zation for both the geometry and the matter fields (see 
e.g. [43! ] for the latter) may even be important in order to 
achieve a satisfactory regularization of any matter field 
divergent contribution. However, along the lines that we 
have commented, the level or even the existence of the 
nonzero energy ground state seems to strongly depend 
on the choice of an evolution parameter and of the par- 
tial observables. This effect may thus restrict possible 
choices within the genuine quantum theory, ruling out 
some selections of an internal time as leading to unphys- 
ical results [44J . 

Another issue which deserves some comments is the 
quantization of the homogeneous degrees of freedom. 
As we have discussed, in the process of quantizing the 
Bianchi I background a la loop, we have applied the 
old quantization prescription for the improved dynam- 
ics, proposed in Refs. [25|,|2g]. Here, the compactness of 
the system protects us against any physical inconsistency 
regarding the dependence on the choice of a fiducial cell 
in the construction. However, it would seem unnatural 
to expect that the compact topology by itself may jus- 
tify using a prescription that is invalid in the noncompact 
cases. Fortunately, the scheme valid for the noncompact 
systems [l3| is sufficiently similar in the aspects that are 
relevant for our analysis, and hence the qualitative re- 
sults reported here must hold in the model built with 
the new scheme as well [23] ]. Nonetheless, the quantita- 
tive results, like the exact bound on the bounce, the scale 
at which the inhomogeneities reach the equilibrium, or 
the critical behavior, are expected to change. Therefore, 
these qualitative results should be viewed as tentative 
only. The main aim of the article is to provide and test 
the methodology, which will be applied in more detailed 
studies of the improved system, exploring and confirming 
the kind of phenomena found here. 

The vacuum Gowdy universes, while simple and use- 
ful to test the formalism and develop the methodology, 
are not fully physical from the viewpoint of observational 
cosmology, since they do not admit near isotropic solu- 
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tions. Furthermore, the degeneracy of the inhomogeneity 
modes qualitatively differs from the observed one, which 
has a completely distinct structure. The first problem 
is solved with the inclusion of matter. For that purpose, 
coupling a massless scalar field to the Gowdy universes is 
a particularly appealing possibility. Building an effective 
theory out of (the already found [3fj|) hybrid quantum de- 
scription of this Gowdy model with matter and repeating 
the dynamical analysis performed here will be the next 
step in our investigations. The second problem, however, 
requires the study of models where the energy level de- 
generacy is that of the spherical harmonics. For that, one 
has to go beyond the nonvacuum Gowdy models, consid- 
ering other more complicated families of spacetimes con- 
taining matter fields and where the inhomogeneities are 
not restricted by the existence of isometries. 
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Appendix A: Different types of data for the 
Monte-Carlo simulations 

In this appendix we explain what kind of universes 
one obtains starting with the different types of initial 
data of Sec. I VII B[ More specifically, we characterize 
the universes with the ratio Vj, /vh , which is small (large) 
for near-homogeneous (inhomogeneities-dominated) uni- 
verses. As can be observed in the first two plots of Fig.[TTl 
when we choose a independent of Vh we get the whole 
range of possible cases, containing both near homoge- 
neous and highly inhomogeneous universes. On the other 

2/3 

hand, when we choose a = av h , our procedure selects 
only a particular region (depending on a) which contains 
just universes peaked around a specific ratio Vb/vh- The 
choice a oc Vh, on the other hand, leads again to solutions 
with a broad range for vj, jvh , but now strongly inhomo- 
geneous universes are generated, owing to the fact that 
a is generically large. Independently of the method of 
generation, the data seem to stay in the same region: a 
well defined patch in the plane A\a m \—(vb/vh)- 
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